/*
   \file  wu.cc
   \brief Implementation of the wu algorithm.

   This file defines the non-inline member functions of the
   gist::wu class. See the header file for more info.
*/

/*
   This file is part of libgist.

   libgist is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2 of the License, or (at your
   option) any later version.

   libgist is distributed in the hope that it will be useful, but WITHOUT
   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
   FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
   for more details.

   You should have received a copy of the GNU General Public License
   along with libgist; if not, write to the Free Software Foundation,
   Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.
*/

/*
   REVISION HISTORY

   $HeadURL: https://libgist.googlecode.com/svn/branches/wu_dev/src/wu.cc $
   $Id: wu.cc 142 2012-02-14 04:02:13Z christian.g.siagian@gmail.com $
*/

//------------------------------ HEADERS --------------------------------

// libgist
#include "wu.hh"

// Boost.Thread
#include <boost/thread/locks.hpp>
#include <boost/thread/mutex.hpp>

// Boost.Timer
#include <boost/timer.hpp>

// Boost.MultiArray
#include <boost/multi_array.hpp>

// Boost.Png_io
#include <boost/gil/extension/io/png_io.hpp>

// Standard C++
#include <iostream>
#include <list>
#include <map>
#include <cstdio>
#include <vector>

// IT++ for Principal Component Analysis
#include <itpp/signal/fastica.h>

//----------------------- NAMESPACE DEFINITION --------------------------

namespace gist {

//----------------------------- CONSTANTS -------------------------------

namespace {

const int CENTRIST_CODEBOOK_METHOD          = 1;
const int SPATIAL_PCA_CENTRIST_METHOD       = 2;   
//const int CURRENT_METHOD                    = CENTRIST_CODEBOOK_METHOD;
const int CURRENT_METHOD                    = SPATIAL_PCA_CENTRIST_METHOD;

// number of reduced feature dimension
const int NUMBER_PCA_COMPONENTS             = 40;
const int NUMBER_CODEBOOK_VOCABULARY        = 200; 

const int NUMBER_SPATIAL_PCA_PYRAMID_LEVELS = 2;
const int NUMBER_CODEBOOK_PYRAMID_LEVELS    = 2;

const int NUMBER_RAW_CENTRIST_VALUES        = 256;

const int L1_NORM_TOTAL                     = 128;

const float SPATIAL_PCA_SCALE_CHANGE        = 0.50;
const float CODEBOOK_SCALE_CHANGE           = 0.75;

const int CENTRIST_WINDOW_SIZE              = 16;
const int CENTRIST_GRID_STEP_SIZE           = 8;

const float UNIFORMITY_THRESHOLD            = 5.0F;

const float VARIANCE_THRESHOLD              = 100.0F;

} // end of local anonymous namespace encapsulating above constants


//-------------------- CENTRIST FEATURES HELPERS ----------------------

namespace {

// CENTRIST (Census Transform) features
// for CODEBOOK: as a codebook histogram
// for spatial PCA: as a PCA output
class centrist_features {

public:

  // Create centrist features
  centrist_features(int num_vals) ;

   // Create centrist features
  centrist_features(std::vector<float> vals) ;

  // get the feature values
  std::vector<float> get_feature_values() ;

  // set a feature value
  void set_value(int index, float val);

  // get the feature size
  uint get_size();

private:
  std::vector<float> its_vals;
} ;


// Constructor: create an array of specified size
centrist_features::centrist_features(int num_vals)
{
  its_vals.clear();
  its_vals.resize(num_vals);
}

// Constructor: create an array of specified size
centrist_features::centrist_features(std::vector<float> vals)
{
  its_vals.clear();
  for(uint i = 0; i < vals.size(); i++)
    its_vals.push_back(vals[i]);
}

// get the feature values
std::vector<float> centrist_features::get_feature_values()
{
  return its_vals;
}

// set feature value at index
void centrist_features::set_value(int index, float val)
{
  assert(index < its_vals.size());
  its_vals[index] = val;
}

// get the vector size
uint centrist_features::get_size()
{
  return its_vals.size();
}

} // end of local anonymous namespace encapsulating above helper

//-------------------- CENTRIST DESCRIPTOR HELPERS ----------------------

// Quick helper to wrap around of a CENTRIST descriptor, which is
// simply a vector of 256 numbers.
namespace {

struct centrist_descriptor {

   static const int SIZE = 256 ;

   // Initialize an empty CENTRIST descriptor, i.e., all zeros.
   centrist_descriptor() ;

   // Return reference to i-th element of the descriptor.
   float& operator[](int i) {return m_values[i] ;}
   const float& operator[](int i) const {return m_values[i] ;}

   // Iterator interface to descriptor values.
   float* begin() {return m_values ;}
   float* end  () {return m_values + SIZE ;}
   const float* begin() const {return m_values ;}
   const float* end  () const {return m_values + SIZE ;}

  // add a value to the vector
  void add(int index, float val = 1.0f);

  // get value
  float get_val(int index);

  // set validity
  void set_validity(bool in_validity);

  // check if the grid is valid (region is distinct) 
  bool is_valid();

  // normalize the histogram
  void normalize();

private :
   float m_values[SIZE] ;
   bool validity;
} ;

// CENTRIST descriptor initialization
centrist_descriptor::centrist_descriptor()
{
   const int n = sizeof m_values / sizeof m_values[0] ;
   std::fill_n(m_values, n, 0) ;

   validity = true;
}

// add value to the value vector
void centrist_descriptor::add(int index, float val)
{
  m_values[index] += val;
}

// get value
float centrist_descriptor::get_val(int index)
{
  return m_values[index];
}

// set validity
void centrist_descriptor::set_validity(bool in_validity)
{
  validity = in_validity;
}

// check validity
bool centrist_descriptor::is_valid()
{
  return validity;
}

// normalize histogram to have 0 mean, 
// remove the first and last entry, and normalize to have unit norm
void centrist_descriptor::normalize()
{
  float sum = 0.0f;
  for(int i = 0; i < SIZE; i++) sum += m_values[i]; sum/= SIZE;

  // normalize histogram to have 0 mean, 
  for(int i = 0; i < SIZE; i++) m_values[i] -= sum;

  // FIXXX:??? remove the first and last entry, 
  //m_values[0] = m_values[255] = 0;

  float sq = 0.0f; 
  for(int i = 0; i < SIZE; i++) sq += m_values[i]*m_values[i]; 

  sq = (sq > VARIANCE_THRESHOLD)?1.0/sqrt(sq):0; 
  for(int i = 0; i < SIZE; i++) m_values[i] *= sq;
}

// compute a centrist descriptor given a window [x1, x2), [y1, y2) 
// compute usually a 16x16 patch 
centrist_descriptor compute_centrist_descriptor
(boost::gil::gray8_image_t ct, int x1,  int x2, int y1, int y2)
{
  boost::gil::gray8_view_t ctv = boost::gil::view(ct);
  boost::gil::gray8_view_t::xy_locator ct_loc = ctv.xy_at(0,0);

  centrist_descriptor cd;

  for(int i = x1; i < x2; i++)
    for(int j = y1; j < y2; j++)
      {
        int val = ct_loc(i,j);
        cd.add(val);
      }

  // normalize the features to zero mean and unit norm
  cd.normalize();

  return cd;
}

} // end of local anonymous namespace encapsulating above helper

//----------------------- CENTRIST DESCRIPTOR GRIDS ------------------

namespace {

// This type facilitates easy iteration over each of the CENTRIST
// descriptors in the grid.
typedef boost::multi_array_ref<const centrist_descriptor, 1> grid_data ;

// The codebook version of centrist (wu, rehg 2011pami) algorithm 
// works by computing centrist descriptors for the
// input image arranged in a regular grid 
// This class is a wrapper around a 2D array to hold them
class centrist_grid {

  // data structure that holds the centrist descriptors for a level 
  // there can be multiple levels depending on the spatial pyramid levels
  boost::multi_array<centrist_descriptor, 2> m_grid ;

  // number of valid/ non-uniform region features
  int m_num_nonuniform_descriptors ;

public:
   // Create a grid of the specified width and height
   centrist_grid(int w, int h) ;

   // Add a centrist descriptor to the specified cell of the grid and update
   // the non-zero descriptors count as appropriate.
   void add(int i, int j, centrist_descriptor d) ;

   // Return the number of non-uniform descriptors in the grid. 
   // This method is meant to be called after all the centrist descriptors for a
   // particular grid have been added.
   int num_nonuniform_descriptors() const {return m_num_nonuniform_descriptors ;}

   // returns an instance of the above type to allow clients
   // to iterate over each of the centrist descriptors contained by this grid
   grid_data data() const {
      return grid_data(m_grid.data(), boost::extents[m_grid.num_elements()]) ;
   }
} ;

// Constructor: create 2D array of specified size
centrist_grid::centrist_grid(int w, int h)
   : m_grid(boost::extents[w][h]), m_num_nonuniform_descriptors(0)
{}

// Add a centrist descriptor to the grid and update non-uniform descriptor count
void centrist_grid::add(int i, int j, centrist_descriptor d)
{
   m_grid[i][j] = d ;
   
   if (d.is_valid())
     ++m_num_nonuniform_descriptors ;
}

// compute the census transform of an image
boost::gil::gray8_image_t compute_census_transform(boost::gil::gray8_image_t img)
{
  boost::gil::gray8_image_t ct(img.dimensions()); 
  boost::gil::gray8_view_t iv = boost::gil::view(img);
  boost::gil::gray8_view_t::xy_locator i_loc = iv.xy_at(0,0);

  boost::gil::gray8_view_t ctv = boost::gil::view(ct);
  boost::gil::gray8_view_t::xy_locator ct_loc = ctv.xy_at(0,0);
    
  const int w = img.width() ;
  const int h = img.height();

  // go through each point 
  for(int y = 1; y < h-1; y++)
    for(int x = 1; x < w-1; x++)
      {
        boost::gil::gray8_pixel_t pxy = i_loc(x,y);

        unsigned char index = 0;
        if(pxy <= i_loc(x-1,y-1)) index += 0x80;
        if(pxy <= i_loc(x-1,y  )) index += 0x40;
        if(pxy <= i_loc(x-1,y+1)) index += 0x20;
        if(pxy <= i_loc(x  ,y-1)) index += 0x10;
        if(pxy <= i_loc(x  ,y+1)) index += 0x08;
        if(pxy <= i_loc(x+1,y-1)) index += 0x04;
        if(pxy <= i_loc(x+1,y  )) index += 0x02;
        if(pxy <= i_loc(x+1,y+1)) index += 0x01;
        ct_loc(x,y) = index;
      }
  return ct;
}

// check whether the region is uniform using 
// intensity standard deviation thresholding
bool check_uniformity
(boost::gil::gray8_image_t img, int x1,  int x2, int y1, int y2)
{
  boost::gil::gray8_view_t img_v = boost::gil::view(img);
  boost::gil::gray8_view_t::xy_locator img_loc = img_v.xy_at(0,0);

  float s1 = 0.0f; float s2 = 0.0f;
  for(int i = x1; i < x2; i++)
    for(int j = y1; j < y2; j++)
      { float val = img_loc(i,j); s1 += val; s2 += val*val; }

  // compute average and variance
  s1 /= ((x2-x1)*(y2-y1)); s2 /= ((x2-x1)*(y2-y1));
  s2 -= s1 * s1;
  
  return (s2 < UNIFORMITY_THRESHOLD);
}

// compute centrist descriptors for 16x16 patches on a grid
// for visual codebook approach
centrist_grid compute_codebook_centrist_grid
(boost::gil::gray8_image_t ct, boost::gil::gray8_image_t img)
{
  const int w = ct.width() ;
  const int h = ct.height();
  
  assert(w >= 3 && h >= 3);
  int half_window = CENTRIST_WINDOW_SIZE/2;
  int wg = (w-half_window-2)/CENTRIST_GRID_STEP_SIZE;
  int hg = (h-half_window-2)/CENTRIST_GRID_STEP_SIZE;
  centrist_grid cg(wg, hg);
  printf("wg: %d hg: %d = %d\n", wg, hg, wg*hg);

  // compute descriptors to each grid location
  for(int i = 0; i < wg; i++)
    {
      for(int j = 0; j < hg; j++)
        {
          int x1 = i*CENTRIST_GRID_STEP_SIZE + 1;
          int x2 = x1 + CENTRIST_WINDOW_SIZE;
          int y1 = j*CENTRIST_GRID_STEP_SIZE + 1;
          int y2 = y1 + CENTRIST_WINDOW_SIZE;
          //printf("x1 %d x2 %d y1 %d y2 %d\n", x1, x2, y1, y2);

          // check for uniformity
          centrist_descriptor cd = 
            compute_centrist_descriptor(ct,x1,x2,y1,y2);          

          bool uniformRegion = check_uniformity(img,x1,x2,y1,y2);
          if(uniformRegion) cd.set_validity(false); // default true

          cg.add(i,j, cd);
        }
      //printf("%d ",i);
    }
  //printf("done compute centrist grid\n");
  return cg;
}

// compute centrist descriptors 
// for a specific spatial pyramid level
centrist_grid compute_spatial_pca_centrist_grid
(boost::gil::gray8_image_t ct, int level)
{
  const int w = ct.width() ;
  const int h = ct.height();
  assert(w >= 3 && h >= 3);

  // for the lowest level
  if(level == 0)
    {
      centrist_grid cg(1, 1);
      cg.add(0,0, compute_centrist_descriptor(ct,0,w-1,0,h-1));
      return cg;
    }

  int num_win_org   = int(pow(2, level));
  int num_win_shift = int(pow(2, level))-1;

  // shift: window/2
  // subtract 2 to take out out of bounds census transform computation
  int w_2 = w-2;
  int h_2 = h-2;
  int win_w_2 = (w_2)/num_win_org/2;
  int win_h_2 = (h_2)/num_win_org/2;

  // just make it 1 dimension
  int num_grid = num_win_org*num_win_org+num_win_shift*num_win_shift;
  centrist_grid cg(num_grid, 1);
  //printf("[%d %d] num grid: %d\n", w,h, num_grid);

  // compute descriptors to each original grid location
  int count = 0;
  for(int i = 0; i < num_win_org; i++)
    {
      for(int j = 0; j < num_win_org; j++)
        {
          int x1 = 1 + int((  i  /float(num_win_org))*w_2);
          int x2 = 1 + int(((i+1)/float(num_win_org))*w_2);
          int y1 = 1 + int(( j   /float(num_win_org))*h_2);
          int y2 = 1 + int(((j+1)/float(num_win_org))*h_2);
          //printf("[%3d] x1 %d x2 %d y1 %d y2 %d\n", count, x1, x2, y1, y2);
          cg.add(count++, 0, compute_centrist_descriptor(ct,x1,x2,y1,y2));
        }
      //printf("%d ",i);
    }

  // also compute the shifted grid location
  // compute descriptors to each original grid location
  int sw = w_2 - win_w_2*2;
  int sh = h_2 - win_h_2*2;
  for(int i = 0; i < num_win_shift; i++)
    {
      for(int j = 0; j < num_win_shift; j++)
        {
          int x1 = 1+win_w_2 + int((  i  /float(num_win_shift))*sw);
          int x2 = 1+win_w_2 + int(((i+1)/float(num_win_shift))*sw);
          int y1 = 1+win_h_2 + int(( j   /float(num_win_shift))*sh);
          int y2 = 1+win_h_2 + int(((j+1)/float(num_win_shift))*sh);
          //printf("[%3d] x1 %d x2 %d y1 %d y2 %d\n", count, x1, x2, y1, y2);
          cg.add(count++, 0, compute_centrist_descriptor(ct,x1,x2,y1,y2));
        }
      //printf("%d ",i);
    }

  //printf("done compute centrist grid\n");
  return cg;
}

// bilinear rescaling
boost::gil::gray8_image_t resize
(boost::gil::gray8_image_t img, float scale)
{
  
  const int orig_w = img.width() ;
  const int orig_h = img.height();

  const int new_w = int(scale*orig_w);
  const int new_h = int(scale*orig_h);
  assert(new_w>0 && new_h>0);

  boost::gil::gray8_image_t result(new_w,new_h);

  // check if same size already
  //if (new_w == orig_w && new_h == orig_h) return img;

  const float sw = float(orig_w) / float(new_w);
  const float sh = float(orig_h) / float(new_h);

  boost::gil::gray8_view_t::iterator dptr = boost::gil::view(result).begin();
  boost::gil::gray8_view_t::iterator sptr = boost::gil::view(img).begin();

  // code inspired from one of the Graphics Gems book:
  /*
    (1) (x,y) are the original coords corresponding to scaled coords (i,j)
    (2) (x0,y0) are the greatest lower bound integral coords from (x,y)
    (3) (x1,y1) are the least upper bound integral coords from (x,y)
    (4) d00, d10, d01, d11 are the values of the original image at the corners
        of the rect (x0,y0),(x1,y1)
    (5) the value in the scaled image is computed from bilinear interpolation
        among d00,d10,d01,d11
  */

  for (int j = 0; j < new_h; ++j)
    {
      const float y = std::max(0.0f, (j+0.5f) * sh - 0.5f);

      const int y0 = int(y);
      const int y1 = std::min(y0 + 1, orig_h - 1);

      const float fy = y - float(y0);

      const int wy0 = orig_w * y0;
      const int wy1 = orig_w * y1;

      for (int i = 0; i < new_w; ++i)
        {
          const float x = std::max(0.0f, (i+0.5f) * sw - 0.5f);

          const int x0 = int(x);
          const int x1 = std::min(x0 + 1, orig_w - 1);

          const float fx = x - float(x0);

          float d00(sptr[x0 + wy0]);
          float d10(sptr[x1 + wy0]);
          float d01(sptr[x0 + wy1]); 
          float d11(sptr[x1 + wy1]);

          float dx0(d00 + (d10 - d00) * fx);
          float dx1(d01 + (d11 - d01) * fx);

          *dptr++ = char( dx0 + (dx1 - dx0) * fy );  // no need to clamp
        }
    }

  return result;
}

// given an input image, 
// compute multiple levels of centrist grid
std::vector<centrist_grid> compute_spatial_pca_centrist_grids
(const boost::gil::gray8_image_t& img)
{
  const int w = img.width() ;
  const int h = img.height();
  std::cout << "CENTRIST image size: " << w << " " << h << "\n";

  // setup how much space is needed given the parameters 
  std::vector<centrist_grid> centrist_grids;

  float scale_change = SPATIAL_PCA_SCALE_CHANGE;
  int curr_pyr_level = NUMBER_SPATIAL_PCA_PYRAMID_LEVELS;

  boost::gil::gray8_image_t uimg = img;

  // FIXXX if use sobel
  //if(USE_SOBEL) { uimg = sobel(img); timg = img; }

  // get the features 
  while(curr_pyr_level >= 0)
    {
      // get the raw census transform features
      boost::gil::gray8_image_t ct = compute_census_transform(uimg);

      // compute centrist grid for this level  
      centrist_grids.push_back
        (compute_spatial_pca_centrist_grid(ct, curr_pyr_level));

      // resize
      uimg = resize(uimg, scale_change);

      // FIXXX if use sobel uimg and timg needs to be rechecked
      //if(USE_SOBEL) { timg = resize(timg, scale_change); uimg = sobel(timg)}

      curr_pyr_level--;
    }

  return centrist_grids ;
}

// given an input image, 
// compute multiple levels of centrist grid
std::vector<centrist_grid> compute_codebook_centrist_grids
(const boost::gil::gray8_image_t& img)
{
  const int w = img.width() ;
  const int h = img.height();
  std::cout << "CENTRIST image size: " << w << " " << h << "\n";

  // setup how much space is needed
  // given the parameters 
  std::vector<centrist_grid> centrist_grids;

  float scale_change = CODEBOOK_SCALE_CHANGE;
  int curr_pyr_level = NUMBER_CODEBOOK_PYRAMID_LEVELS;

  boost::gil::gray8_image_t uimg = img;

  // FIXXX if use sobel
  //if(USE_SOBEL) { uimg = sobel(img); timg = img; }

  // get the features 
  while(curr_pyr_level >= 0)
    {
      // get the raw census transform features
      boost::gil::gray8_image_t ct = compute_census_transform(uimg);

      // compute centrist grid for this level  
      centrist_grids.push_back
        (compute_codebook_centrist_grid(ct, uimg));

      // resize
      uimg = resize(uimg, scale_change);

      // FIXXX if use sobel timg and uimg needs to be rechecked
      //if(USE_SOBEL) { timg = resize(timg, scale_change); uimg = sobel(timg)}

      curr_pyr_level--;
    }

  return centrist_grids ;
}

} // end of local anonymous namespace encapsulating above helpers

//---------CENTRIST GRIDS AND FEATURES LISTS ---------------------------

namespace {

// This data structure holds all the centrist grids for a 
// particular image category. It is used in training mode.
typedef std::list<std::vector<centrist_grid> > centrist_grids_list ;

// Quick helper to create a new list of centrist grids containing
// the grid passed in to the function.
inline centrist_grids_list new_list(const std::vector<centrist_grid>& c)
{
   centrist_grids_list list ;
   list.push_back(c) ;
   return list ;
}

// This data structure holds all the centrist features for a 
// particular image category. It is used in training mode.
typedef std::list<centrist_features> centrist_features_list ;

// Quick helper to create a new list of centrist features containing
// the features passed in to the function.
inline centrist_features_list new_list(const centrist_features& c)
{
   centrist_features_list list ;
   list.push_back(c) ;
   return list ;
}

} // end of local anonymous namespace encapsulating above helpers


//---------------------- WU IMPLEMENTATION ------------------------

// The wu class hides its entire implementation behind an opaque
// object. This section defines this implementation class.
namespace {

// In training mode, we will accumulate the centrist grids for the
// training images in a map indexed by the category labels.
typedef std::map<std::string, centrist_grids_list> centrist_grids_map ;

// In training mode, we then produce centrist features for the
// training images after further processing the centrist grids; 
// the features are stored in a map indexed by the category labels.
typedef std::map<std::string, centrist_features_list> centrist_features_map ;

// This class actually implements the Wu algorithm. The interface
// class defined in the header file, for the most part, simply forwards
// its API to this one.
class wu_impl {

  // stores the centrist grids to store raw centrist features
  centrist_grids_map cg_map;

  // This data structure holds the final centrist features for the input
  // images arranged according to their categories.
  centrist_features_map c_map ;
  
  // As we allow multithreaded clients, we will have to synchronize
  // access to the above members.
  mutable boost::mutex m_mutex ;

public:
  // Constructor
  wu_impl() ;
  
  // Add a new centrist grids given its category
  void add_centrist_grids
  (const std::string& cat, const std::vector<centrist_grid>&) ;
  
  // Add a new centrist features given its category
  void add_centrist_features
  (const std::string& cat, const centrist_features&) ;

  // compute pca training 
  void compute_pca_training();

  // compute k-means codebook training
  void compute_codebook_training();

  // train place classifier
  void train_place_classifier();

} ;

// Constructor
wu_impl::wu_impl()
{}

// lock map of centrist grids and add new features given its category
void wu_impl::add_centrist_grids
(const std::string& cat, const std::vector<centrist_grid>& cgrids)
{
  boost::lock_guard<boost::mutex> lock(m_mutex) ;
  centrist_grids_map::iterator m = cg_map.find(cat) ;
  if (m == cg_map.end())
    cg_map.insert(std::make_pair(cat, new_list(cgrids))) ;
  else
    m->second.push_back(cgrids) ;
}

// lock map of centrist features and add new features given its category
void wu_impl::add_centrist_features
(const std::string& cat, const centrist_features& cf)
{
  boost::lock_guard<boost::mutex> lock(m_mutex) ;
  centrist_features_map::iterator m = c_map.find(cat) ;
  if (m == c_map.end())
    c_map.insert(std::make_pair(cat, new_list(cf))) ;
  else
    m->second.push_back(cf) ;
}

// 
void wu_impl::compute_pca_training()
{
  boost::timer t ;   
  std::cout << "WU: PCA and ICA dimension reduction\n";

  // calculate the number of descriptors in training
  uint num_examples = 0;
  uint num_labels = cg_map.size();  
  centrist_grids_map::iterator cgptr = cg_map.begin();  
  for (uint l = 0; l < num_labels; ++l)
    {            
      centrist_grids_list cglist = (*cgptr++).second;

      uint list_size = cglist.size(); 
      centrist_grids_list::iterator cglptr = cglist.begin();
      for(uint i = 0; i < list_size; i++)
        {
          std::vector<centrist_grid> cg = (*cglptr++);
          uint num_grid = cg.size();
          for(uint j = 0; j < num_grid; j++)
            {
              uint num_descriptor =  cg[j].num_nonuniform_descriptors();

              //std::cout << l << " " << i << " " << j << " " << num_descriptor << "\n"; 
              num_examples += num_descriptor;
            }
        }
    }
  printf("num_examples: %d\n", num_examples);

  // fill in the signal to be fed to the PCA/ICA engine
  itpp::mat  signal(centrist_descriptor::SIZE,num_examples);
  
  cgptr = cg_map.begin(); uint curr_cd = 0;
  for (uint l = 0; l < num_labels; ++l) 
    {      
      centrist_grids_list cglist = (*cgptr++).second;

      uint list_size = cglist.size(); 
      centrist_grids_list::iterator cglptr = cglist.begin();
      for(uint i = 0; i < list_size; i++)
        {
          std::vector<centrist_grid> cg_vec = (*cglptr++);
          uint num_grid = cg_vec.size();
          for(uint j = 0; j < num_grid; j++)
            {
              uint num_descriptor =  cg_vec[j].data().size();
              grid_data cg = cg_vec[j].data();
              grid_data::iterator cg_ptr = cg.begin();
              for(uint k = 0; k < num_descriptor; k++)
                {
                  centrist_descriptor cd = *cg_ptr;
                  if(!cd.is_valid()) continue;

                  //std::cout << l << " " << i << " " << j << " " << k << "\n"; 
                  for (int c = 0; c < centrist_descriptor::SIZE; ++c) 
                    {
                      signal(c,curr_cd) = cd.get_val(c);
                      //if(float(cd.get_val(c)) > 0.0f)
                      //  printf("[%3d]: %f ", c, float(cd.get_val(c)));
                    }
                  cg_ptr++;
                  curr_cd++;
                  //std::cout << "please enter a character\n"; 
                  //char str[256]; std::cin >> str;
                }
            }
        }
    }

  double s2 = t.elapsed();
  std::cout << "moving data " << s2 << " seconds\n" ;

  std::cout << "mat: " << signal.rows() << " " << signal.cols() << " " 
            << NUMBER_PCA_COMPONENTS << "\n";  

  // compute ICA signals
  itpp::Fast_ICA fastICA(signal);
  fastICA.set_nrof_independent_components(NUMBER_PCA_COMPONENTS);
  // fastICA.set_non_linearity(FICA_NONLIN_TANH); // or FICA_NONLIN_POW3 
  fastICA.set_approach( FICA_APPROACH_DEFL ); // Much faster
  fastICA.separate();

  // the separating matrix
  itpp::mat W = fastICA.get_separating_matrix();
  std::cout << "Save the separating map W: " << W.rows() 
            << " " << W.cols() << "\n";
  uint n_ic = fastICA.get_nrof_independent_components();
  std::cout << "num components: " << n_ic << "\n";

  // store the separated signals
  // same as: ica_signal = W*signal; 
  itpp::mat ica_signal = fastICA.get_independent_components();
  std::cout << "ica_signal: " << ica_signal.rows() << " " << ica_signal.cols() << "\n";

  cgptr = cg_map.begin(); curr_cd = 0;
  for (uint l = 0; l < num_labels; ++l) 
    {
      // get the label and the grid list
      std::string label = (*cgptr).first;
      centrist_grids_list cglist = (*cgptr++).second;

      uint list_size = cglist.size(); 
      centrist_grids_list::iterator cglptr = cglist.begin();
      for(uint i = 0; i < list_size; i++)
        {
          std::vector<float> temp_vals;

          std::vector<centrist_grid> cg_vec = (*cglptr++);
          uint num_grid = cg_vec.size();
          for(uint j = 0; j < num_grid; j++)
            {
              uint num_descriptor =  cg_vec[j].data().size();
              grid_data cg = cg_vec[j].data();
              grid_data::iterator cg_ptr = cg.begin();
              for(uint k = 0; k < num_descriptor; k++)
                {
                  centrist_descriptor cd = *cg_ptr;
                  if(cd.is_valid())
                    {
                      //std::cout << l << " " << i << " " << j << " " << k << "\n"; 
                      // bool same = true;
                      // itpp::mat  signal2(centrist_descriptor::SIZE,1);
                      // for (int c = 0; c < centrist_descriptor::SIZE; ++c) 
                      //   {
                      //     signal2(c,0) = cd.get_val(c);
                      //   }
                      // itpp::mat ica_signal2 = W*signal2; 
                      
                      for (uint c = 0; c < n_ic; ++c) 
                        {
                          temp_vals.push_back(ica_signal(c,curr_cd));
                          //signal(c,curr_cd) = cd.get_val(c);

                          // if(ica_signal(c, curr_cd) != ica_signal2(c,0))
                          //   same = false;
                          //printf("[%3d]: %f = %f \n", 
                          //       c, ica_signal(c, curr_cd), ica_signal2(c,0));
                        }

                      // if(!same)
                      //   {
                      //     std::cout << "please enter a character\n"; 
                      //     char str[256]; std::cin >> str;
                      //   }
                    }
                  else
                    {
                      // insert zero vector
                      for (uint c = 0; c < n_ic; ++c) 
                        temp_vals.push_back(0.0f);
                    }
                  cg_ptr++;
                  curr_cd++;
                }
            }

          // create a centrist_features
          //printf("size: %d", int(temp_vals.size()));
          add_centrist_features(label, centrist_features(temp_vals)) ;
        }
    }

  // to run it: 
  double s3 = t.elapsed();
  std::cout << "PCA/ICA takes: " << s3 << " seconds\n" ;
}

// 
void wu_impl::compute_codebook_training()
{
  boost::timer t ;   
  std::cout << "WU: create vocabulary codebook\n";

}

//
void wu_impl::train_place_classifier()
{
  // train 
  std::cout << "\nWU: train SVM place classifier FIXXX: TODO\n";

  // char buffer[200];
  // sprintf(buffer,"image_%d.png", curr_pyr_level);
  // boost::gil::png_write_view(buffer, boost::gil::view(timg));

  //std::cout << "please enter a character\n"; 
  //char str[256]; std::cin >> str;

  // after PCA, etc
  //float curr_weight  = 1.0f;
  //curr_weight *= (1.0f/scale_change/scale_change);

  // gist_features_map gf_map = pimpl->get_gist_features_map() ;  
}

} // end of local anonymous namespace encapsulating above implementation details

//------------------------ WU INTERFACE ---------------------------

// Initialization
wu::wu()
   : p_impl(new wu_impl)
{}

// collect Training data
void wu::
data(const std::string& cat, const boost::gil::gray8_image_t& img)
{
  wu_impl* pimpl = reinterpret_cast<wu_impl*>(p_impl) ;
  
   // compute and store the centrist features
  if(CURRENT_METHOD == SPATIAL_PCA_CENTRIST_METHOD)
    pimpl->add_centrist_grids(cat, compute_spatial_pca_centrist_grids(img)) ;

  else if(CURRENT_METHOD == CENTRIST_CODEBOOK_METHOD)
    pimpl->add_centrist_grids(cat, compute_codebook_centrist_grids(img)) ;
}

// training data
void wu::train()
{
  wu_impl* pimpl = reinterpret_cast<wu_impl*>(p_impl) ;
  
  // perform PCA on each grid region histogram
  if(CURRENT_METHOD == SPATIAL_PCA_CENTRIST_METHOD)
    pimpl->compute_pca_training();

  // or perform codebook training on 16x16 centrist patch
  else if(CURRENT_METHOD == CENTRIST_CODEBOOK_METHOD)
    pimpl->compute_codebook_training();
  
  // train place classifiers
  pimpl->train_place_classifier();

  // FIXXX: Lazebnik does local helper copy "class copy"
  
  //std::cout << "please enter a character\n"; 
  //char str[256]; std::cin >> str;
}

// Clean-up
wu::~wu()
{
   delete reinterpret_cast<wu_impl*>(p_impl) ;
}

//-----------------------------------------------------------------------

} // end of namespace encapsulating this file's definitions

/* So things look consistent in everyone's emacs... */
/* Local Variables: */
/* indent-tabs-mode: nil */
/* End: */
